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The fluctuations in wind power entering an electrical grid (Irish grid) were analyzed and found 
to exhibit correlated fluctuations with a self-similar structure, a signature of large-scale correlations 
in atmospheric turbulence. The statistical structure of temporal correlations for fluctuations in 
generated and forecast time series was used to quantify two types of forecast error: a timescale error 
(e T ) that quantifies the deviations between the high frequency components of the forecast and the 
generated time series, and a scaling error (e^) that quantifies the degree to which the models fail to 
predict temporal correlations in the fluctuations of the generated power. With no a priori knowledge 
of the forecast models, we suggest a simple memory kernel that reduces both the timescale error 
(e T ) and the scaling error (e^). 


INTRODUCTION 

Renewable power generation, unlike conventional 
power, exhibits variability owing to natural fluctuations 
in the energy source [Tj. Wind power, in particular, 
shares spectral features of the turbulent wind from which 
it derives energy m ■ This variability in power output 
adds a cost to renewable power 00 that is absent in con¬ 
ventional power sources. Whereas distributed wind farms 
are expected to smooth the fluctuations 0 , power enter¬ 
ing the electrical grid still exhibits large amplitude fluctu¬ 
ations |Tj. Large ramps in power fluctuations present the 
possibility of grid destabilization [5] and blackout, a con¬ 
stant source of concern for system operators 0 (9j. This 
risk increases the cost of operating reserves 110! needed 
on standby to return a grid back to operation in the event 
of failure. Naturally, forecast models constitute essential 
tools in estimating the magnitude of fluctuations before¬ 
hand and in planning for the optimal operating reserves 
required on call. Yet, no standards for forecast accuracy 
currently exist El- 

Extant works on wind power forecast error, ranging 
from the turbine to the grid scale, focus on modeling the 
forecast error distribution [7214171 . Since a probability 
distribution is time-independent, it contains no informa¬ 
tion on temporal error variations. Several studies have 
considered the dependence of the mean and the variance 
of the error on the duration for which the power is pre¬ 
dicted (ranging from minutes to hours) [£ TTl . Other 
works have considered the different distributions of er¬ 
rors for mean power over different durations mmm- 
However, none of these studies account for the fluctuation 
correlations of atmospheric turbulence [21] transferred to 
the generated power in the analysis of forecast error nor 
for the temporal correlations of the errors. 

Whereas power fluctuations at the scale of an individ¬ 
ual turbine 0 and a wind farm 00 have been shown 


to exhibit self-similar scaling, such fluctuations from in¬ 
dividual wind farms are expected to smooth out before 
they enter the electrical grid. Using data from the Irish 
grid operator EIRGRID i.22], we show that wind power 
entering the grid exhibits correlated fluctuations with a 
self-similar structure. Such scaling points to large-scale 
correlations in atmospheric turbulence influencing the ag¬ 
gregate wind power entering the grid. 

In this article, we exploit these correlations at the grid 
scale and draw upon the Statistical Theory of Hydrody¬ 
namic Turbulence to quantify two types of forecast er¬ 
ror. The first is a timescale error (e T ) that quantifies the 
timescales over which the forecast models fail to predict 
high-frequency power fluctuations. This timescale error 
sets a bound on the numerical resolution of forecast mod¬ 
els and would already be known to system operators who 
own and run the forecast models. However, details of the 
models are not available to potential customers in energy 
spot markets [23] who could use this information to fac¬ 
tor in the risk associated with a non-supply of promised 
wind power by producers. The second type of error we 
quantify is a scaling error (ej) that establishes a differ¬ 
ence in the self-similar scaling of fluctuations as observed 
for actual generated power vis a vis the power that was 
forecast to be generated. This error could be potentially 
useful to model developers, and if such an error results 
from large-scale correlations in atmospheric turbulence, 
incorporating them into models is not subject to limi¬ 
tations arising from numerical resolution. Having estab¬ 
lished the errors, we then employ a simple memory kernel 
upon the forecast time series and show that the errors can 
be easily reduced with a minimal computational cost. 

Two raw time series are provided by EIRGRID: the 
wind power generated nationwide across all Ireland en¬ 
tering the grid p g {t), and the power forecast by EIR- 
GRID’s modelsp/(t) for the same period. The time series 
sampled at 15 minute intervals span a five-year period 
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FIG. 1: (color online) a) Raw time series (for 10 days) of 
the generated power p g (t) (black open circles), forecast power 
Pf{t) (red open squares), and the instantaneous forecast error 
Pd(t) (blue open triangles) in megawatts (MW). Every third 
data point is plotted for easy visibility, b) The probability 
density function of the raw instantaneous forecast error II (pd) 
(solid black circles) has exponentially decaying, fat tails rela¬ 
tive to a Gaussian distribution (solid black line) of the same 
mean and standard deviation as II(p t j). 

(2009 — 2014). No other information that permits mean¬ 
ingful data decomposition is available; forecast models 
employed, the number of wind farms feeding the grid, 
their location, date of commission or date of scheduled 
and unscheduled outages, etc. are all unknown. Despite 
the lack of further information, any quantifiable trends 
revealed are of immediate use to system operators in es¬ 
timating operating reserves and in accounting for fluc¬ 
tuations that could potentially destabilize the grid 123] , 
Furthermore, such information [25j is potentially useful 
to customers in energy spot markets [25] . 

Raw time series for the generated p g (t) and forecast 
power pf(t), and their instantaneous difference Pd{t) = 


Pf(t) — Pg{t), which we define as the instantaneous fore¬ 
cast error, are shown in Fig. Eh for a 10-day period, 
permitting a few immediate qualitative observations. 
Firstly, p g (t) exhibits correlated fluctuations. Secondly, 
Pf(t ), while closely following p g (t), misses the high fre¬ 
quency (relative to the sampling rate of time series) com¬ 
ponents. The instantaneous forecast error pd(t) exhibits 
correlated fluctuations with a coefficient of variation- 
standard deviation/mean = 116.24/22.9 ~ 5, implying 
large magnitude fluctuations in pd{t) (i.e., a broad dis¬ 
tribution) . 

DISTRIBUTION OF FORECAST ERROR 

As a point of comparison with prior works EOT , we 
note that the probability density function (PDF) of the 
raw instantaneous forecast error Il(pd) (fig- EE) exhibits 
fat exponential tails that decay slower than a Gaussian 
function of the same mean and standard deviation as 
n(p rf ). Indeed, a PDF with exponential tails may be 
expected for reasons detailed in the following. Being a 
scalar product of an instantaneous force (/(f)) and ve¬ 
locity the statistics of temporal variation in power 

are determined by the product of two random variables 
p(t) = v(t) ■ f(t) [25]. The statistics of the product 
(Z = XY ) of two normally distributed random variables 
(X and Y) was first studied by C. Craig [27] (henceforth 
referred to as Craig’s-XY distribution). An asymptotic 
analysis reveals that Craig’s-XY distribution is logarith¬ 
mically singular about zero with exponentially decaying 
tails, and its asymmetry (skewness) depends on the in¬ 
stantaneous cross-correlation between X and Y (28,, [29]. 
Whereas this asymptotic analysis is possible only when 
X and Y are Gaussian, the structure of Craig’s XY- 
distribution itself is more generally observed, even when 
X and/or Y are non-Gaussian [28]. 

The basic structure of Craig’s-XY distribution is ex¬ 
pected for the PDF of power p(t) (e.g., compare Fig. 2c in 
[3] with Fig. 2 in [29]) and its estimation error Sp(t). Sup¬ 
pose for a single wind turbine, the errors in estimating or 
forecasting velocity and force are Sv(t) and 6f(t), respec¬ 
tively, then p(t) + 5p(t ) = (v(t) + Sv(t)) ■ (/(f) + 6f(t)). 
Expanding the RHS permits decomposition into power 
p{t) = v(t) ■ f(t) and its error 

Sp(t) = v(t) ■ 6f(t) + 6v(t) ■ f(t) + Sv(t) ■ Sf(t). (1) 

Consequently, 5p(t) is a random variable whose statis¬ 
tics are determined by the sum of three terms (Eq. EJ , 
each being a product of two random variables. One ex¬ 
pects 8p(t) to exhibit features of Craig’s-XY distribu¬ 
tion irrespective of whether or not v, 5v, /, and 5f are 
Gaussian. In fact, given that the velocity distribution of 
atmospheric turbulence is known to follow the Weibull 
distribution ESI IS], a numerical approach may become 
necessary. 
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The power forecast error statistics can be easily scaled 
from the turbine to the grid scale. If M wind farms feed 
power to the grid, and the ith farm has Ni turbines, the 
cumulative error in estimating wind power then equals 
the instantaneous forecast error p d (t) for the grid, shown 
in Fig. UK and is given by: 

p d (t) = x?i i xyi 1 6p j (t). ( 2 ) 

Summing the power error statistics over all turbines 
(across all farms) causes an averaging of fluctuations, 
starting with the most probable ones that occur around 
zero, thus smoothing the logarithmic singularity. All 
these generic features are readily observed for II (p d ) in 
Fig. [T], Beyond contributing to the extant literature 
UMII, the structure of II(pd) provides no useful infor¬ 
mation for our analysis and will not be discussed fur¬ 
ther. Understanding temporal variability (fluctuations) 
and uncertainty (error) requires analysis of the temporal 
evolution of the distributions, their moments and multi¬ 
point temporal correlation functions. We therefore pro¬ 
ceed through a statistical analysis of the temporal cor¬ 
relations in the fluctuating time series for generated and 
forecast power. 


DATA ANALYSIS 


The time series was analyzed in two stages, with trends 
in the series being identified in the first stage, followed 
by an analysis of the fluctuations around the trends in 
the second stage. Trend removal permits a focus on sys¬ 
tematic differences between p g (t) and Pf(t), ignoring dif¬ 
ferences due to new wind farms and seasonal variability 
of the wind power. Trend identification was performed 
such that the cross-correlation between the generated and 
forecast power trends was maximal. We used the fast 
Fourier transform (FFT) for each of the series and defined 
the trends by inverting the transform using only the fre¬ 
quencies with maximal amplitudes. The number of max¬ 
imal amplitudes was set by the requirement of the high¬ 
est cross-correlation between p g (t) and Pf{t). Keeping 
the zero frequency (to preserve the signal mean) and five 
more frequencies resulted in a peak cross-correlation of 
0.9904 between the generated and forecast power trends 
(Fig#)- These respective trends were subtracted from 
the raw time series. We denote the de-trended gener¬ 
ated power by Pc(t), forecast power by Pp(f) and their 
instantaneous difference by Pn{t) = Pf{P) — Pg(P)- 
The characteristic fluctuation timescales for the de¬ 
trended time series were first computed from their re¬ 
spective autocorrelation functions defined as: 


C X ( T ) = (-Py(t) ~ p x)(Px(t + t) - P x ) 

(Px(t)-p^r 


(3) 


where Px is a time-average subtracted from the signal 
(de-trending does not render a zero signal mean since the 
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FIG. 2: (color online) a) The five-year trend for p g (t) (black 
solid line) and P;(t) (red dashed line) is subtracted from the 
raw time series in subsequent analysis, b) Log-linear scale: 
autocorrelation functions Cg{t) (open black circles), Cf(t) 
(open red squares) and Cd(t) (open blue triangles) for Pg(I), 
Pf (t) and Po(t), respectively, exhibit exponential decorrela¬ 
tion with respective characteristic timescales obtained from 
fit to data of tg = 80.94 points (20.24 hours), tf = 81 points 
(20.24 hours) and td = 25.86 points (~ 6.5 hours). Every 
third data point is plotted for easy visibility. 


zero frequency component was preserved). The subscript 
X should be replaced with G for generated power, F for 
forecast power, and D for instantaneous forecast error, 
respectively. The three autocorrelation functions (fig. §>) 
exhibit exponential decay for short times with a data 
fit following the functional form C'jy(r) ~ Axe~^ T / Tx \ 
where Ax — 1.0, owing to Cx(t) being normalized, 
and t x represents the characteristic decorrelation time 
for each time series, yielding tq = 80.94 data points 
(~ 20.24 hours) for generated power, tf = 81 points 
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(also ~ 20.24 hours) for forecast power, and td = 25.86 
points (~ 6.5 hours) for instantaneous forecast error. The 
de-trended series were also split into independent time se¬ 
ries of shorter duration (l/8th of the original temporal 
duration). Autocorrelation functions computed for these 
windowed data did not reveal a measurable difference 
in the characteristic decay time tx\ deviations were ap¬ 
parent only for long-term behavior spanning a week (or 
longer timescales) when the decorrelation had already 
occurred. The correlation time of high frequency fluctu¬ 
ations (< 20 hours) is much shorter than the slow varying 
trend (over months to years). Hence the de-trending pro¬ 
tocol (in particular, the number of maximal amplitudes) 
does not influence the analysis to follow-a fact verified 
and reported upon later. 

Autocorrelation functions for the generated (Cg(t)) 
and forecast (Cf(t)) power exhibit nearly identical scal¬ 
ing and the same characteristic decay timescales (tq = 
t f = 20.24 hours), suggesting the accurate capture of 
correlations in generated power by the forecast models. 
Yet, the autocorrelation function Cfo(r) for instanta¬ 
neous forecast error Pd{P) informs us that some corre¬ 
lations are not captured. In particular, we qualitatively 
know that P F (t) misses the high frequency components 
of and they end up in Po(t), thereby contribut¬ 

ing to its two-point correlator. This correlation deficit 
suggests that higher order moments of the two-point cor¬ 
relator are necessary to capture the statistical structure 
of the missing fluctuations. 

TEMPORAL STRUCTURE FUNCTIONS 

Statistical analysis of higher order correlations is a 
well-developed, mature tool within the Statistical The¬ 
ory of Hydrodynamic Turbulence in which higher or¬ 
der two-point correlators are studied through Structure 
Functions. Kolmogorov’s theory of 1941 (K41) [32] lays 
the foundation for structure functions through the cel¬ 
ebrated “4/5th law”: S^r) = ((Aw|| (r)) 3 ) = ((v\\(R + 
r) — U||(-R)) 3 ) = — |er, where the third moment of lon¬ 
gitudinal velocity differences (((At>||(r)) 3 )) between two 
points spatially separated by a longitudinal distance r, 
is proportional to the product of the average turbulent 
dissipation rate (e) and the longitudinal spacing r [551 . 

The nth order structure function encodes all cross¬ 
terms up to order n of the two-point correlator for a given 
stationary signal. The physical relevance of structure 
functions may be appreciated by considering a stationary, 
fluctuating signal x(t) with zero mean. The difference 
between the two values of this signal taken time r apart 
(A x(t) = x{t + r) — x(t)) is collected at various windows 
(of duration r) along the time series. A x(t) is therefore 
a random variable with statistics of its own, and the nth 
order structure function defined as S n (r) = ((Ax(r))") 
is the nth moment for its PDF n(Aa;(r)). The moment 



FIG. 3: (color online) Structure functions of order n = 1 — 10 
(red solid circles) and their power-law fits (black solid lines) 
for (a) generated power S (r) and (b) forecast power (t) 
plotted versus r in log-log scale exhibit self-similar scaling 
Sn (r) oc "A" (X is G for generated and F for forecast power). 
The scaling is robust for (a) the generated power over 1.4 
decades (40 time steps), (b) In contrast, for forecast power, 
the first- and second-order structure functions exhibit scaling 
up to r = 40 time steps, but for n > 2, no scaling is observed 
for r < 10 time steps. Self-similar scaling is restored over a 
limited range of timescales 10 < r < 40. 


S n (r ) varies with the time difference r between signals, 
and its scaling (if any) reveals temporal variations in the 
statistical structure of fluctuations in the signal to the 
nth order. 

Tails of the PDF n(Ax(r)) exert themselves with the 
increasing order n of the structure function, thus necessi¬ 
tating more data to resolve higher order structure func¬ 
tions. A weak test for resolving the nth order structure 
function involves splitting the time series into smaller 
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windows and testing for identical scaling on the trun¬ 
cated series. However, this test only assures stationarity 
of the statistics. A strong test for the ability to resolve 
the nth order structure function requires that first, the 
moment’s integrand (Ax) n n(Ax) —► 0 as |Ax| — ► oo [53] 
(required due to finiteness of data), and second, the PDF 
n(Ax) should decay faster than 1/| Ax | n+1 for | Ax| oo 
or else the integral / (Ax) ra n(Ax) dx would diverge for 
large |Ax| [35] (test for existence of a PDF’s nth mo¬ 
ment). Whereas the two conditions are not indepen¬ 
dent, the second condition is theoretical and does not 
depend upon the available statistics. When conducting 
data analysis, even when the second condition is satis¬ 
fied, insufficient data can lead to noise and prevent the 
integrand (Ax) n n(Ax) from satisfactorily converging to 
zero. The first condition is therefore dependent on finite¬ 
ness of data. Based on both weak and strong tests, we 
conclude that the EIRGRID data can resolve structure 
functions up to order n = 12; however, only results up to 
n = 10 are presented. 

Since even-order structure functions take only positive 
values, they converge faster than ones with odd order. To 
overcome this distinction between odd and even orders, 
we compute the nth order structure function of the abso¬ 
lute value of differences: S* (r) = (| Px(t + r) — Px{t)\ n ) 
where subtraction of mean Px(t + r) and Px(t) is as¬ 
sumed. While ensuring the same convergence rate for 
even- and odd-order statistics, it also collates all data 
in the positive quadrant permitting easy visualization. 
Analysis of fractional-order structure functions allows 
better testing for anomalous scaling [36:. Fractional- 
order structure functions are only defined for absolute 
values of signal differences [35]-another reason why we 
calculate structure functions of absolute differences. We 
also calculated the structure functions of orders n = 
0.1 — 0.9 in steps of 0.1. 

RESULTS 

Figure [3] plots the structure functions of order n = 
1 — 10 (fractional-order structure functions are calcu¬ 
lated but not shown) for the absolute value of signal 
differences of the generated power |A(Pg(t))| (fig. 
and forecast power |A(Pp(r))| (Fig. |3 Jd). Self-similar 
or power-law scaling is observed for the generated power 
structure functions over 1.4 decades spanning r < 40. 
Scaling over the same temporal range is also observed 
for the forecast power structure functions of order n = 1 
and 2. For n > 2, no scaling is observed for timescales 
r < 10. The scaling is restored over a limited range of 
timescales 10 < r < 40 (0.4 decades in time). 

Self-similar scaling of the temporal structure functions 
implies a relationship of the form: 

Sn( T ) A nT C * (4) 


where (,* is the scaling exponent. For simple mono¬ 
fractal scaling, (£ oc n. However, fluctuations with a 
multi-fractal character exhibit a nonlinear dependence of 
the scaling exponent (,* with respect to n. Super- (sub-) 
linear variation of C,* versus n implies temporal expan¬ 
sion (compression) of fluctuations [33 ■ Scaling exponents 
for all the structure functions were computed from the 
log derivative, Cn = d ^loglr)' 1 ' 1 > which provides a more 
reliable estimate of the exponent than a power-law fit 
[351 |33] . The pre-factor A* in Eq. [ 4 ] is subsequently 
obtained from fit to data. In Fig. [3j all the data (red 
solid circles) were divided by A% such that all fits (solid 
black lines) commence from both mantissa (t) and ordi¬ 
nate (jS^(t)) at unity, for an easy comparison of (n with 
order n. All the data in Fig. [3] Fig.[4]a,, and Fig. 0 3 there¬ 
fore follow the scaling relation: S*(t) oc (A* = 1). 


The scaling in Fig.[3]reveals higher-order temporal cor¬ 
relations at work in the EIRGRID data. The absence of 
scaling for (t) for n > 2 at timescales r < 10 con¬ 
firms the qualitative observation made in Fig. [T]i that 
forecast models do not capture high frequency fluctu¬ 
ations. More importantly, Fig. HJ 3 ascribes a precise 
bound on the time (t = 10, 2.5 hours) out to which the 
high frequency fluctuations are missed. Finally, scaling 
presence for S^(t), n = 1,2 explains the close agree¬ 
ment between the autocorrelation functions Cg(t) and 
Cf(t) and their identical characteristic decay times, tq 
and Tp, observed in Fig. i 3 - This is to be expected 
on the grounds that the second-order structure function 
Sz{t) = ((Ax(t )) 2 ) = (x(f+T) 2 )+(x(f) 2 )—2(x(f)x(f+r)) 
shares a direct correspondence with the autocorrelation 
function where the cross-term is identical to the numer¬ 
ator of Eq. [3| The failure of (t) for n > 2 to capture 
high frequency fluctuations out to r = 10 reveals one type 
of forecast error in the models; we call this the timescale 
error e T . 


Before proceeding to the second type of error aris¬ 
ing from scaling mismatch, we define the cross-structure 
function X£ g (t) = (\P F (t + T) - P G (t)\ n ). X^ g (t) rep¬ 
resents nth order moments for the PDF of the relative 
magnitude of fluctuations between Pq (t) and P F (t + r ), 
and their cross-terms correspond to higher-order two- 
point cross-correlators between the generated and fore¬ 
cast power. This function is plotted in Fig. [4]i. Again, 
we notice that scaling is absent at early times (r < 10), 
and restored at later times (10 < r < 40). We note 
that X£ g (t) exhibits no scaling for n = 1 and 2, un¬ 
like the forecast structure functions (Fig. [3jo). Although 
Sn ( r ) exhibits scaling for order n = 1 and 2, its expo¬ 
nent (n'i this scaling deficit is reflected in X^ g (t) 

for n = 1 and 2. 
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FIG. 4: (color online) a) Log-log scale: cross-structure func¬ 
tions Xn G (j) versus r (red solid circles) exhibit no scaling at 
early times r < 10, with scaling restored for 10 < r < 40. 
Black solid lines are power-law fits to data within the scaling 
regime, b) Scaling exponent versus the order of struc¬ 
ture function n for generated G (red solid circles), forecast F 
(blue solid squares) and modified forecast M (black solid tri¬ 
angles) structure functions, and cross-structure functions FG 
(green solid inverted triangles), and their respective second- 
order polynomial fits: solid red line for small dashed blue 
line for , medium dashed black line for and long dashed 
green line for C,n°■ 


DISCUSSION 

Having established the various structure functions, we 
now consider the behavior of their scaling exponents Cn 
(X = G for generated, F for forecast and FG for the 
cross-structure function). Figure [ijr plots versus 
the order n together with their polynomial fits to the 
quadratic order. ( G = 10 -2 4- 0.67n — 0.013?r 2 scales 
almost linearly (mono-fractal) with a small, but mea¬ 


surable, quadratic deviation towards multi-fractal behav¬ 
ior. The exponent = 0.007 + 0.8n — 0.025n 2 exhibits 
a slightly more pronounced quadratic deviation (multi 
fractal behavior) relative to ( G . On the other hand, 
(^ G = 10~ 2 +0.54n — 0.006n 2 scales almost linearly with 
n, implying mono-fractal scaling. 

We now consider the measurement error for the afore¬ 
mentioned scalings. Firstly, given that all de-trending 
protocols suffer from an ad hoc choice of a de-trending 
timescale, we tested the scalings for dependence on the 
de-trending procedure by varying the number of max¬ 
imal amplitudes. Ignoring the condition for maximal 
cross-correlation between p g (t) and Pf(t), the number of 
maximal amplitudes contributing to the trends was var¬ 
ied. The scalings were invariant up to the inclusion of 15 
maximal amplitudes into the trend, beyond which, co¬ 
efficients for the polynomial fits started varying in the 
second decimal place. Having ascertained the robustness 
of our choice for the five maximal amplitudes at which the 
cross-correlation peaks, we focused on a second source of 
scaling measurement error, namely statistical variability. 
Since the scalings are analyzed up to r = 100 data points, 
the de-trended time series were split into eight indepen¬ 
dent windows (each with 21912 data points), and the 
structure functions were re-computed for each window. 
The variation in the log derivative (C,f = d ) 

for the eight independent measurements was taken as the 
possible scatter in the scaling estimation, thereby pro¬ 
viding a confidence interval for the polynomial fits. The 
scatter was found to be ± 0.01 in both the measured 
value of (£ and the corresponding polynomial fits (for 
each of the polynomial coefficients) for each of the eight 
independent datasets, revealing that the polynomial fits 
were meaningful only to the linear order for Q G and Cn G ■ 
The quadratic-order polynomial coefficient for de¬ 
spite being larger than the scatter of ±0.01, is not useful 
owing to the fact that the corresponding quadratic terms 
for ( G and ^ G are smaller than the scatter magnitude. 

Despite qualitatively observing a quadratic deviation 
for in Fig. [4 }d, our inability to ascribe significance to it 
arises from the fact that the multi-fractal component (de¬ 
viation from linear scaling) of the scalings is miniscule. 
This is significant in light of several studies that have 
demonstrated multi-fractal scaling for wind power fluc¬ 
tuations at the turbine 001 and farm scales 02 ■ Turbu¬ 
lence theory traces the source of multi-fractal behavior to 
intermittent fluctuations that can arise from two sources 
in the atmospheric context. The first, known as inter¬ 
nal intermittency, occurs at the small scales of turbulent 
flow. These intermittent fluctuations would be naturally 
reflected in the power generated at the turbine and farm 
scales. However, when adding together power generated 
by geographically distant wind farms, internal intermit¬ 
tency should smooth out 0D] since it is a small-scale ef¬ 
fect and cannot extend across geographically distributed 
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wind farms. Furthermore, the sampling interval (15 min¬ 
utes) for EIRGRID data is not expected to resolve any 
effects that may arise from internal intermittency, which 
occur at much shorter timescales (high frequencies). 

The second source of intermittency, known as external 
intermittency, occurs at the edge of any free-stream [JT| 
and arises in the atmospheric context due to coupling be¬ 
tween the atmospheric boundary layer turbulence and a 
co-moving weather system ]21j . External intermittency, 
which can be experienced in the form of wind gusts, is of 
greater relevance in the present analysis as it can both 
correlate distributed farms through the weather system 
and occur at timescales longer than the 15-minute sam¬ 
pling interval for EIRGRID data. The nearly fractal scal¬ 
ing of informs us that both internal and external in¬ 
termittency are being smoothed to the point of rendering 
grid-level power fluctuations almost mono-fractal. 

The self-similar scaling of S^(t) over several hours 
does strongly point to the influence of large-scale turbu¬ 
lent structures on power fluctuations at the grid level. 
The 20-hour characteristic decorrelation time (tq) for 
generated power in Fig. n, if taken as the large eddy 
turnover time of atmospheric turbulence, also lends cre¬ 
dence to such an argument. Finally, independent proof 
in support of this argument also comes from Katzenstein 
et al. 20 ] who show that an individual wind farm ex¬ 
hibits f ~ 5 / 3 (/ being the frequency) scaling for the wind 
power spectrum (equivalent to r 2 / 3 scaling of the second- 
order structure function in the time domain). However, 
as wind power from various farms is summed, the spec¬ 
trum steepens (please see Fig. 3 in [40]). Such spectral 
steepening can be clearly attributed to the smoothing 
of high frequency (short timescale) fluctuations corre¬ 
sponding to small eddies. But the low frequency (long 
timescale) fluctuations corresponding to large-scale ed¬ 
dies lose no power spectral density, clearly indicating the 
influence of large-scale turbulent structures. 

We finally consider the forecast error due to the scaling 
mismatch. We define the scaling error as ~ Cn ■ 

Under this definition, if the time series for forecast and 
generated power were identical, then S^(t) = 
implying , and therefore = 0. Another typical 

case arises if forecast models fail completely, resulting in 
a flat time series with no fluctuations, ^ = 0 resulting 
in an error er = ~Cn ■ Using the polynomial fits for 
(see Fig. Eh) to linear order, we obtain = (7 x 
10~ 2 + 0.8n) - (HT 2 + 0.67n) = -0.003 + 0.13n. This 
can be cross-validated against the difference G = 

(10 _2 +0.67n) —(10 _2 + 0.54?r) = 0.13n. Since —>■ 0 as 

n —> 0, the 0th order term falling within the scatter may 
be taken to be zero. Both estimates of error are identical 
in linear order (e^ = 0.013n). 

The analysis thus far demonstrates the importance of 
temporal correlations in wind power and their role in es¬ 
timating forecast errors. It is reasonable to ask whether 
this knowledge could help improve the forecast time se- 



FIG. 5: (color online) a) Log-linear scale: 7 op t versus order of 
structure function shows no improvement for n < 4 but shows 
better agreement for n > 4 with an abrupt change observed 
in 7o P t at n = 4. b) Log-log scale: Structure functions S?? (t) 
versus r (red solid circles) for the modified forecast time series 
show considerable improvement over their counterparts S„ (r) 
in Fig. [ 3 J 3 . 

ries, despite having no knowledge of the models em¬ 
ployed. In particular, to capture the short-term corre¬ 
lations missed by the forecast, we introduce a modified 
forecast that is based on the original forecast, convoluted 
with an exponentially decaying memory kernel derived 
from the generated power time series. The modified fore¬ 
cast power is given by P]\i(t) = f 0 PF(r)e~' r ^ t ~ T ^dT. 

The memory duration (I/ 7 ) was chosen so as to min¬ 
imize the relative difference between the structure func¬ 
tions of the generated and forecast power. As expected 
(as shown earlier, the low-order structure functions of the 
generated and forecast power are very similar), we found 
that the optimal 7 varies with the order of the structure 









function. For n < 4, the memory-modified forecast shows 
no improvement in the agreement between and - 
For n > 4, the modified forecast exhibits better agree¬ 
ment with the structure functions of the generated power 
as shown in Fig. 5ji. The optimal 7 (7opt) was found to 
be 74 « 1.06 and 710 « 0.37, as shown in Fig. plot¬ 
ted in log-linear scale to show the variation in 7 opt for 
n > 4. The simple scheme, suggested here, not only tries 
to rectify the timescale error e T , but also attempts to 
statistically align the temporal correlations by improv¬ 
ing the scaling error e^. 

As is apparent from Fig. [5 Jd, the structure functions 
(Sn(r) = {|APm(t)I”)) for modified forecast time series 
are substantially improved over their unmodified coun¬ 
terpart (fig. §)). Firstly, scalings are restored at high 
frequencies (r < 10), thus rendering the timescale er¬ 
ror irrelevant. More importantly, the scaling itself is 
improved as is evident from Fig. |4|), revealing ^ = 
0.01 + 0.7n — 0.007n 2 . To linear order, the scaling er¬ 
ror — (n= 0.7 n — 0.67n = 0.03n, a considerable 

improvement over the original forecast time series. Be¬ 
ing computationally inexpensive, and given that spinning 
and non-spinning reserves must act within 10 minutes 
of failure, with replacement reserves acting within 20-60 
minutes, there are tangible benefits to incorporating such 
a memory kernel into models to monitor instabilities in 
real-time. Furthermore, it might be possible to improve 
the forecast models using different parameterizations of 
the regional climate models or weather models. 

SUMMARY 

In summary, wind power exhibits significant tempo¬ 
ral correlations even at the grid level, where fluctuations 
are expected to average out [7] as power is fed from 
geographically distributed wind farms. Previous stud¬ 
ies show that the temporal correlations of the wind are 
essential to studying wind-generated large-scale ocean 
currents |42j : a similar appreciation of large-scale cor¬ 
relations in atmospheric turbulence within the context 
of wind power is called for. Fluctuations, albeit pos¬ 
ing a problem to system operators, possess a statistical 
structure through temporal correlations, which could be 
exploited to quantitatively analyze the error in forecast 
models. The technique proposed here is only limited by 
the sampling rate of the time series. Beyond potentially 
serving as a standard for quantifying wind-power forecast 
accuracy, it could have applications for any renewable en¬ 
ergy source with temporally correlated fluctuations pos¬ 
sessing a statistical structure. 
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